options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # not see parameters str1..., str2,... using regex as exc='[str1,str2,...]' 
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(4,1),mar=c(1,5,1,1))
      drw[,1,] |> plot(type='l',xlab='',ylab=pr)
      drw[,2,] |> plot(type='l',xlab='',ylab=pr)
      drw[,3,] |> plot(type='l',xlab='',ylab=pr)
      drw[,4,] |> plot(type='l',xlab='',ylab=pr)
      par(mar=c(3,5,3,3))
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}
fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)

mcmc=goMCMC(mdl,data,smp,wrm)

mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')


y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')


m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')
}

explanatory variable obs.x have noise

xN(x0.sx),yN(b0+b1*x0,s)

ex8-0.stan

normal regression without explanatory variable’s noise

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}

ex9.stan

normal regression with explanatory variable’s noise

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x10;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
  real<lower=0> sx;
  vector[N] x0;
}  
model{
  for(i in 1:N){
    x[i]~normal(x0[i],sx);
    y[i]~normal(b0+b1*x0[i],s);
  }
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x0[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] x1;
  vector[N1] y1;
  for(i in 1:N1){
    x1[i]=normal_rng(x10[i],sx);
    m1[i]=b0+b1*x10[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x0=runif(n,0,20)
x=rnorm(n,x0,2)
y=rnorm(n,x0*2+5,2)
qplot(x,y)

n1=10

#explanatory variable do not has noise
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan') 
fn(mdl,data)
## Initial log joint probability = -33279.9 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       25      -29.4414    0.00414791   0.000567954           1           1       53    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
##  variable estimate
##    lp__     -29.44
##    b0         6.69
##    b1         1.80
##    s          2.64
##    m0[1]     47.65
##    m0[2]      8.21
##    m0[3]     26.38
##    m0[4]     27.20
##    m0[5]     29.29
##    m0[6]     22.88
##    m0[7]     21.40
##    m0[8]     41.24
##    m0[9]     18.12
##    m0[10]    16.17
##    m0[11]    30.09
##    m0[12]    20.15
##    m0[13]    28.04
##    m0[14]    36.42
##    m0[15]    42.89
##    m0[16]    42.80
##    m0[17]    11.65
##    m0[18]     8.08
##    m0[19]     9.86
##    m0[20]    10.35
##    y0[1]     43.65
##    y0[2]      4.97
##    y0[3]     26.89
##    y0[4]     26.51
##    y0[5]     32.87
##    y0[6]     24.22
##    y0[7]     23.61
##    y0[8]     47.38
##    y0[9]     17.91
##    y0[10]    15.97
##    y0[11]    30.39
##    y0[12]    18.75
##    y0[13]    29.75
##    y0[14]    36.71
##    y0[15]    44.16
##    y0[16]    46.78
##    y0[17]    16.28
##    y0[18]     8.91
##    y0[19]     4.77
##    y0[20]     6.08
##    m1[1]      8.08
##    m1[2]     12.47
##    m1[3]     16.87
##    m1[4]     21.27
##    m1[5]     25.66
##    m1[6]     30.06
##    m1[7]     34.46
##    m1[8]     38.86
##    m1[9]     43.25
##    m1[10]    47.65
##    y1[1]      5.09
##    y1[2]      8.37
##    y1[3]     15.28
##    y1[4]     20.41
##    y1[5]     24.67
##    y1[6]     28.37
##    y1[7]     34.87
##    y1[8]     43.52
##    y1[9]     44.01
##    y1[10]    47.75
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -30.13 -29.75 1.40 1.12 -32.99 -28.65 1.01      604      841
##    b0       6.63   6.59 1.22 1.20   4.70   8.63 1.01      266      392
##    b1       1.81   1.80 0.10 0.09   1.64   1.97 1.01      360      680
##    s        3.01   2.95 0.58 0.50   2.23   4.10 1.00      891      986
##    m0[1]   47.74  47.72 1.47 1.39  45.40  50.18 1.00     1059     1286
##    m0[2]    8.16   8.11 1.15 1.12   6.35  10.06 1.01      268      391
##    m0[3]   26.39  26.38 0.72 0.66  25.22  27.59 1.00     1061      875
##    m0[4]   27.21  27.20 0.73 0.67  26.04  28.42 1.00     1294      838
##    m0[5]   29.31  29.30 0.76 0.72  28.11  30.56 1.00     1835     1063
##    m0[6]   22.88  22.84 0.72 0.69  21.76  24.09 1.01      556      730
##    m0[7]   21.39  21.37 0.73 0.71  20.25  22.62 1.01      455      642
##    m0[8]   41.31  41.30 1.17 1.08  39.41  43.24 1.00     1587     1208
##    m0[9]   18.10  18.07 0.80 0.76  16.85  19.44 1.01      348      465
##    m0[10]  16.15  16.13 0.85 0.80  14.80  17.61 1.01      316      503
##    m0[11]  30.11  30.10 0.78 0.74  28.88  31.38 1.00     1996      954
##    m0[12]  20.14  20.12 0.75 0.72  18.95  21.37 1.01      403      563
##    m0[13]  28.06  28.05 0.74 0.68  26.87  29.27 1.00     1523      961
##    m0[14]  36.47  36.47 0.97 0.92  34.88  38.05 1.00     2127     1240
##    m0[15]  42.96  42.95 1.24 1.16  40.95  45.03 1.00     1414     1180
##    m0[16]  42.87  42.86 1.24 1.15  40.87  44.93 1.00     1424     1164
##    m0[17]  11.61  11.57 1.01 0.95  10.02  13.29 1.01      280      408
##    m0[18]   8.02   7.98 1.16 1.13   6.20   9.93 1.01      268      391
##    m0[19]   9.81   9.76 1.08 1.03   8.14  11.60 1.01      273      391
##    m0[20]  10.31  10.26 1.06 1.01   8.65  12.07 1.01      274      395
##    y0[1]   47.80  47.80 3.43 3.27  42.23  53.30 1.00     1614     1877
##    y0[2]    8.15   8.10 3.27 3.06   2.96  13.87 1.00     1126     1624
##    y0[3]   26.35  26.37 3.17 2.98  21.21  31.41 1.00     2127     2100
##    y0[4]   27.19  27.22 3.15 2.90  22.07  32.56 1.00     1861     1963
##    y0[5]   29.35  29.29 3.15 2.97  24.28  34.60 1.00     1769     1764
##    y0[6]   22.86  22.89 3.14 3.05  17.69  27.78 1.00     1907     1891
##    y0[7]   21.36  21.35 3.26 3.00  16.24  26.75 1.00     1782     1896
##    y0[8]   41.31  41.31 3.35 3.04  35.87  46.69 1.00     1913     1971
##    y0[9]   18.17  18.12 3.21 3.06  12.99  23.38 1.00     1789     1971
##    y0[10]  16.30  16.37 3.16 3.12  11.05  21.33 1.00     1611     1698
##    y0[11]  30.30  30.26 3.19 3.06  25.09  35.53 1.00     1956     1749
##    y0[12]  20.11  20.10 3.13 3.00  14.90  25.10 1.00     1680     1956
##    y0[13]  28.14  28.21 3.16 2.94  22.87  33.31 1.00     1728     1941
##    y0[14]  36.57  36.56 3.26 3.17  31.30  41.83 1.00     1937     1972
##    y0[15]  42.99  42.96 3.24 3.13  37.73  48.38 1.00     1778     1821
##    y0[16]  42.90  42.95 3.27 3.15  37.41  48.19 1.00     2021     1700
##    y0[17]  11.61  11.67 3.23 3.16   6.12  16.82 1.00     1653     1681
##    y0[18]   8.10   7.97 3.28 3.06   2.74  13.50 1.00     1352     1564
##    y0[19]   9.75   9.82 3.16 2.94   4.43  14.81 1.00     1345     1681
##    y0[20]  10.30  10.19 3.24 3.06   5.14  15.78 1.01      973     1672
##    m1[1]    8.02   7.98 1.16 1.13   6.20   9.93 1.01      268      391
##    m1[2]   12.44  12.39 0.98 0.92  10.89  14.08 1.01      285      423
##    m1[3]   16.85  16.82 0.83 0.79  15.54  18.26 1.01      325      493
##    m1[4]   21.26  21.24 0.73 0.70  20.12  22.49 1.01      449      615
##    m1[5]   25.68  25.66 0.72 0.66  24.53  26.89 1.00      912      829
##    m1[6]   30.09  30.08 0.78 0.74  28.85  31.36 1.00     1991      954
##    m1[7]   34.50  34.50 0.90 0.86  33.04  35.97 1.00     2270     1207
##    m1[8]   38.91  38.92 1.07 0.99  37.16  40.67 1.00     1862     1190
##    m1[9]   43.33  43.31 1.26 1.18  41.29  45.43 1.00     1379     1210
##    m1[10]  47.74  47.72 1.47 1.39  45.40  50.18 1.00     1059     1286
##    y1[1]    8.15   8.15 3.34 3.18   2.75  13.56 1.00      947     1287
##    y1[2]   12.38  12.31 3.11 3.04   7.37  17.49 1.00     1367     1644
##    y1[3]   16.78  16.77 3.11 3.03  11.71  22.06 1.00     1503     1785
##    y1[4]   21.14  21.13 3.20 3.01  15.89  26.19 1.00     2127     1872
##    y1[5]   25.77  25.86 3.11 3.00  20.75  30.78 1.00     1802     1847
##    y1[6]   30.21  30.20 3.15 3.02  25.10  35.38 1.00     2052     1919
##    y1[7]   34.61  34.65 3.22 3.14  29.25  39.92 1.00     1685     1692
##    y1[8]   38.86  38.86 3.30 3.15  33.49  44.09 1.00     1967     1895
##    y1[9]   43.36  43.33 3.31 3.15  37.97  49.02 1.00     2010     1946
##    y1[10]  47.72  47.69 3.35 3.23  42.43  53.44 1.00     1611     1497

#explanatory variable has noise
x10=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x10=x10)

mdl=cmdstan_model('./ex9.stan')
mcmc=goMCMC(mdl,data,wrm=500,smp=1000)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 1 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 2 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 2 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 3 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 4 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 3 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 4 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 1 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 2 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 3 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 1 finished in 0.4 seconds.
## Chain 2 finished in 0.4 seconds.
## Chain 3 finished in 0.4 seconds.
## Chain 4 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 4 finished in 0.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.4 seconds.
## Total execution time: 0.6 seconds.
seeMCMC(mcmc,'[m,x]')
##  variable   mean median    sd   mad     q5   q95 rhat ess_bulk ess_tail
##    lp__   -27.18 -29.35 10.87 10.30 -42.65 -9.01 1.04       71      107
##    b0       7.07   6.88  1.56  1.31   4.80 10.37 1.11       23       24
##    b1       1.77   1.78  0.13  0.11   1.50  1.96 1.12       23        8
##    s        1.91   1.96  1.02  1.21   0.33  3.51 1.11       26       10
##    sx       1.17   1.19  0.60  0.65   0.20  2.12 1.09       31       93
##    x0[1]   21.89  22.04  1.01  1.08  20.10 23.26 1.02      157      603
##    x0[2]   -0.15   0.11  1.25  1.05  -3.12  1.39 1.11       24       12
##    x0[3]   10.15  10.18  0.84  1.02   8.93 11.45 1.05       67     1545
##    x0[4]   11.51  11.48  0.69  0.53  10.39 12.71 1.03     1714     1673
##    x0[5]   12.54  12.54  0.66  0.54  11.41 13.63 1.01     1087     1902
##    x0[6]    9.17   9.11  0.69  0.61   8.13 10.36 1.02      145     1927
##    x0[7]    6.31   6.40  1.53  1.86   3.50  8.36 1.11       25        9
##    x0[8]   19.38  19.29  0.82  0.69  18.07 20.84 1.07       34       88
##    x0[9]    6.10   6.15  0.77  0.68   4.83  7.34 1.06       47      275
##    x0[10]   4.02   4.22  1.24  1.30   1.55  5.63 1.11       25        9
##    x0[11]  12.32  12.32  0.78  0.89  11.12 13.53 1.02      175     2177
##    x0[12]   8.20   8.22  0.80  0.84   6.96  9.50 1.03      306     1433
##    x0[13]  11.60  11.61  0.68  0.60  10.52 12.71 1.01      514     2123
##    x0[14]  17.58  17.47  1.08  1.16  16.08 19.61 1.10       27       34
##    x0[15]  21.00  20.79  1.13  1.00  19.55 23.58 1.11       24       83
##    x0[16]  21.33  21.11  1.30  1.33  19.67 23.94 1.12       22       35
##    x0[17]   3.33   3.27  0.84  0.93   2.08  4.71 1.05       64      338
##    x0[18]   0.93   0.94  0.89  0.78  -0.66  2.34 1.10       27       31
##    x0[19]   1.91   1.94  0.86  0.73   0.25  3.29 1.08       30       13
##    x0[20]   2.94   2.87  0.96  1.07   1.61  4.53 1.03      138     1395
##    m0[1]   45.74  45.48  1.93  2.08  43.18 49.18 1.03      149     1014
##    m0[2]    6.91   6.68  1.60  1.56   4.68  9.64 1.02      204     1656
##    m0[3]   25.05  24.90  1.46  1.54  22.90 27.47 1.02      179     1805
##    m0[4]   27.43  27.48  1.26  1.05  25.32 29.49 1.04     3948     1821
##    m0[5]   29.23  29.23  1.22  1.02  27.19 31.27 1.02     3530     1961
##    m0[6]   23.31  23.35  1.24  1.02  21.19 25.27 1.03      795     1743
##    m0[7]   18.33  18.15  2.38  3.20  15.08 22.07 1.04       80      861
##    m0[8]   41.29  41.28  1.35  1.14  39.08 43.51 1.02     2963     2674
##    m0[9]   17.90  17.91  1.25  1.01  15.91 19.99 1.04     2510     1933
##    m0[10]  14.27  14.09  1.76  2.12  11.87 17.13 1.03      131      840
##    m0[11]  28.87  28.74  1.40  1.43  26.86 31.28 1.02      204     2092
##    m0[12]  21.58  21.70  1.49  1.69  19.13 23.71 1.03      144     2182
##    m0[13]  27.59  27.48  1.23  1.08  25.63 29.69 1.04     1142     1670
##    m0[14]  38.10  38.28  1.66  1.80  35.16 40.34 1.03      135      899
##    m0[15]  44.12  44.35  1.59  1.51  41.27 46.32 1.02      211      896
##    m0[16]  44.69  44.98  1.85  1.91  41.36 47.18 1.02      170      789
##    m0[17]  13.00  13.15  1.51  1.60  10.47 15.11 1.02      171     1659
##    m0[18]   8.78   8.94  1.42  1.21   6.27 10.88 1.01      679     1943
##    m0[19]  10.51  10.66  1.36  1.19   8.24 12.64 1.03      813     1773
##    m0[20]  12.31  12.54  1.78  1.95   9.31 14.78 1.03      139     1147
##    y0[1]   45.76  45.12  2.89  2.35  41.99 51.17 1.01      460     1591
##    y0[2]    6.95   6.47  2.70  2.09   3.04 11.95 1.01      767     2474
##    y0[3]   25.05  24.63  2.52  2.03  21.39 29.63 1.02      753     1704
##    y0[4]   27.46  27.48  2.49  1.85  23.39 31.65 1.04     4033     2651
##    y0[5]   29.21  29.19  2.51  1.90  25.13 33.44 1.04     3788     1906
##    y0[6]   23.36  23.44  2.51  1.92  19.12 27.32 1.04     3393     2606
##    y0[7]   18.33  17.62  3.26  3.10  14.32 24.31 1.03      135     1810
##    y0[8]   41.35  41.27  2.56  1.97  37.14 45.72 1.03     3763     3186
##    y0[9]   17.86  17.83  2.57  1.87  13.71 22.34 1.03     3113     2782
##    y0[10]  14.30  13.69  2.82  2.34  10.59 19.48 1.01      364     1498
##    y0[11]  28.87  28.42  2.62  2.00  25.03 33.81 1.02      928     2099
##    y0[12]  21.59  22.03  2.62  2.14  16.87 25.26 1.02      679     2011
##    y0[13]  27.56  27.43  2.46  1.85  23.48 31.88 1.04     3814     3358
##    y0[14]  38.18  38.70  2.70  2.21  33.16 41.93 1.01      416     2176
##    y0[15]  44.12  44.56  2.66  2.10  39.33 47.93 1.02      897     2533
##    y0[16]  44.71  45.33  2.85  2.23  39.38 48.62 1.01      433     1892
##    y0[17]  12.98  13.39  2.67  2.13   8.19 16.89 1.01      593     2684
##    y0[18]   8.71   9.01  2.60  1.99   4.05 12.67 1.02     2333     2817
##    y0[19]  10.54  10.74  2.58  1.99   5.97 14.66 1.05     3089     2202
##    y0[20]  12.28  12.90  2.80  2.30   7.05 15.97 1.02      386     2314
##    m1[1]    8.44   8.25  1.47  1.23   6.27 11.51 1.11       24       24
##    m1[2]   12.76  12.62  1.19  1.03  10.95 15.16 1.11       24       26
##    m1[3]   17.08  16.98  0.94  0.85  15.62 18.84 1.10       26       33
##    m1[4]   21.40  21.36  0.75  0.75  20.20 22.57 1.07       39      299
##    m1[5]   25.73  25.75  0.66  0.63  24.61 26.73 1.02      336     2311
##    m1[6]   30.05  30.00  0.71  0.66  28.90 31.20 1.01      664     2036
##    m1[7]   34.37  34.37  0.89  0.90  33.03 35.80 1.04       63     1123
##    m1[8]   38.70  38.73  1.13  1.08  36.90 40.49 1.07       35      125
##    m1[9]   43.02  43.07  1.40  1.29  40.68 45.18 1.09       29       50
##    m1[10]  47.34  47.39  1.68  1.52  44.38 49.91 1.10       27       39
##    x1[1]    0.74   0.75  1.33  0.98  -1.48  2.97 1.02     4338     1139
##    x1[2]    3.24   3.23  1.31  0.98   1.05  5.45 1.02     4268     1673
##    x1[3]    5.65   5.65  1.31  0.97   3.47  7.81 1.02     4201     1602
##    x1[4]    8.10   8.10  1.31  0.96   5.99 10.26 1.02     4137     1727
##    x1[5]   10.52  10.54  1.31  0.97   8.23 12.71 1.02     4032     1117
##    x1[6]   13.00  13.00  1.32  0.95  10.81 15.23 1.02     3733     1580
##    x1[7]   15.43  15.43  1.29  0.93  13.22 17.55 1.02     4008     1294
##    x1[8]   17.87  17.88  1.28  0.97  15.77 19.93 1.02     3976     2072
##    x1[9]   20.31  20.32  1.28  0.92  18.16 22.42 1.02     3946     1624
##    x1[10]  22.77  22.78  1.34  0.96  20.48 24.98 1.02     3972     1511
##    y1[1]    8.41   8.41  2.68  2.40   4.01 12.28 1.05       55     2207
##    y1[2]   12.85  12.84  2.45  2.31   8.86 16.42 1.04       88     2845
##    y1[3]   17.08  17.14  2.39  2.11  13.17 20.80 1.03      115     2504
##    y1[4]   21.40  21.49  2.33  1.81  17.48 25.14 1.01     1359     2477
##    y1[5]   25.71  25.82  2.28  1.69  21.78 29.47 1.03     2863     2433
##    y1[6]   30.11  30.04  2.22  1.64  26.46 33.86 1.04     3095     2120
##    y1[7]   34.38  34.29  2.35  1.85  30.64 38.51 1.01     1188     2880
##    y1[8]   38.71  38.65  2.39  2.10  35.02 42.74 1.02      149     2805
##    y1[9]   42.98  42.98  2.56  2.36  39.12 47.33 1.03      102     2725
##    y1[10]  47.37  47.38  2.75  2.54  43.23 51.84 1.04       70     2183

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
x1=mcmc$draws('x1')
smx1=summary(x1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=smx1$median,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

outlier

ex10.stan

objective variable have outlier, y~cauchy(b0+b1*x,s)

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~cauchy(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=cauchy_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=cauchy_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,9)
y=rnorm(n,x*2+5,1)
x[1]=3
y[1]=25
qplot(x,y)

n1=10

x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan') 
fn(mdl,data)
## Initial log joint probability = -29125.9 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       28      -33.6454   0.000773142   0.000329433           1           1       59    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -33.65
##    b0         6.46
##    b1         1.85
##    s          3.26
##    m0[1]     12.01
##    m0[2]     14.93
##    m0[3]     15.25
##    m0[4]     21.16
##    m0[5]     16.82
##    m0[6]     18.91
##    m0[7]     14.07
##    m0[8]     18.98
##    m0[9]     10.35
##    m0[10]    13.62
##    m0[11]    11.20
##    m0[12]    13.32
##    m0[13]    23.08
##    m0[14]     8.65
##    m0[15]    22.81
##    m0[16]    11.18
##    m0[17]    22.49
##    m0[18]    18.66
##    m0[19]    17.31
##    m0[20]    15.28
##    y0[1]      6.71
##    y0[2]     16.07
##    y0[3]     17.33
##    y0[4]     27.00
##    y0[5]     15.21
##    y0[6]     19.38
##    y0[7]     12.65
##    y0[8]     20.67
##    y0[9]      7.39
##    y0[10]     8.21
##    y0[11]     8.83
##    y0[12]    16.14
##    y0[13]    24.05
##    y0[14]     8.64
##    y0[15]    22.00
##    y0[16]    12.15
##    y0[17]    23.97
##    y0[18]    18.96
##    y0[19]    16.16
##    y0[20]    16.52
##    m1[1]      8.65
##    m1[2]     10.25
##    m1[3]     11.86
##    m1[4]     13.46
##    m1[5]     15.06
##    m1[6]     16.67
##    m1[7]     18.27
##    m1[8]     19.87
##    m1[9]     21.48
##    m1[10]    23.08
##    y1[1]      7.33
##    y1[2]     11.79
##    y1[3]     11.44
##    y1[4]     17.09
##    y1[5]     15.45
##    y1[6]     18.33
##    y1[7]     17.94
##    y1[8]     16.61
##    y1[9]     17.38
##    y1[10]    18.59
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -34.04 -33.73 1.22 1.14 -36.46 -32.63 1.00      614     1031
##    b0       6.71   6.64 2.01 1.92   3.27  10.14 1.02      303      337
##    b1       1.81   1.81 0.36 0.35   1.21   2.41 1.01      358      422
##    s        3.70   3.61 0.66 0.60   2.77   4.95 1.00     1238     1161
##    m0[1]   12.14  12.14 1.12 1.08  10.26  14.02 1.01      344      596
##    m0[2]   14.99  14.98 0.85 0.84  13.64  16.45 1.00      735     1353
##    m0[3]   15.31  15.29 0.84 0.82  13.98  16.74 1.00      848     1353
##    m0[4]   21.09  21.10 1.32 1.30  18.94  23.22 1.00      981     1270
##    m0[5]   16.84  16.83 0.85 0.84  15.49  18.29 1.00     2109     1587
##    m0[6]   18.89  18.89 1.01 1.01  17.23  20.58 1.00     2030     1492
##    m0[7]   14.15  14.13 0.90 0.85  12.73  15.68 1.00      512      994
##    m0[8]   18.96  18.96 1.02 1.02  17.28  20.67 1.00     1962     1492
##    m0[9]   10.51  10.49 1.36 1.29   8.23  12.83 1.01      313      424
##    m0[10]  13.71  13.69 0.94 0.89  12.19  15.30 1.01      447      783
##    m0[11]  11.34  11.31 1.23 1.18   9.30  13.43 1.01      325      473
##    m0[12]  13.42  13.41 0.97 0.91  11.84  15.06 1.01      416      730
##    m0[13]  22.97  22.98 1.62 1.60  20.29  25.54 1.00      736     1035
##    m0[14]   8.85   8.82 1.63 1.55   6.10  11.66 1.01      305      339
##    m0[15]  22.71  22.72 1.58 1.55  20.10  25.21 1.00      757     1071
##    m0[16]  11.32  11.29 1.24 1.18   9.27  13.42 1.01      324      473
##    m0[17]  22.39  22.40 1.53 1.50  19.87  24.82 1.00      789     1095
##    m0[18]  18.64  18.64 0.99 0.98  17.03  20.29 1.00     2163     1537
##    m0[19]  17.33  17.31 0.87 0.87  15.92  18.80 1.00     2254     1569
##    m0[20]  15.34  15.32 0.84 0.83  14.01  16.77 1.00      859     1352
##    y0[1]   12.19  12.23 4.01 3.89   5.57  18.86 1.00     1559     1789
##    y0[2]   15.03  15.11 3.85 3.66   8.82  21.25 1.00     1882     2015
##    y0[3]   15.30  15.25 3.79 3.51   9.06  21.68 1.00     2017     2014
##    y0[4]   21.26  21.34 3.88 3.83  14.86  27.34 1.00     1984     1933
##    y0[5]   16.63  16.67 3.84 3.61  10.31  22.73 1.00     2036     1786
##    y0[6]   18.92  18.96 3.97 3.81  12.38  25.22 1.00     1923     1862
##    y0[7]   14.04  13.99 3.94 3.72   7.58  20.55 1.00     1773     1845
##    y0[8]   18.92  19.00 3.99 3.78  12.08  25.46 1.00     1890     1911
##    y0[9]   10.47  10.50 3.98 3.81   3.94  17.22 1.01     1168     1856
##    y0[10]  13.70  13.57 3.89 3.86   7.40  20.02 1.00     1986     2014
##    y0[11]  11.40  11.29 3.90 3.88   4.85  17.77 1.00     1223     1452
##    y0[12]  13.38  13.29 3.88 3.83   7.34  19.65 1.00     1499     1812
##    y0[13]  22.92  22.82 3.99 3.96  16.53  29.80 1.00     1832     1835
##    y0[14]   8.84   8.80 4.07 3.88   2.19  15.52 1.00     1023     1971
##    y0[15]  22.78  22.86 4.11 3.93  16.09  29.33 1.00     1621     1881
##    y0[16]  11.27  11.29 4.00 3.86   4.80  18.09 1.00     1320     1693
##    y0[17]  22.34  22.36 4.12 4.14  15.70  28.96 1.00     1556     1830
##    y0[18]  18.72  18.72 3.86 3.67  12.51  25.17 1.00     2232     1933
##    y0[19]  17.43  17.42 3.81 3.72  11.19  23.56 1.00     1900     1802
##    y0[20]  15.36  15.33 3.88 3.71   9.37  21.64 1.00     1930     1808
##    m1[1]    8.85   8.82 1.63 1.55   6.10  11.66 1.01      305      339
##    m1[2]   10.42  10.40 1.37 1.30   8.11  12.77 1.01      312      420
##    m1[3]   11.98  11.98 1.14 1.09  10.06  13.91 1.01      340      544
##    m1[4]   13.55  13.54 0.95 0.89  12.01  15.17 1.01      430      744
##    m1[5]   15.12  15.11 0.84 0.84  13.79  16.56 1.00      779     1316
##    m1[6]   16.69  16.68 0.84 0.82  15.36  18.14 1.00     1926     1594
##    m1[7]   18.26  18.26 0.95 0.93  16.73  19.85 1.00     2250     1629
##    m1[8]   19.83  19.83 1.13 1.12  17.95  21.70 1.00     1377     1369
##    m1[9]   21.40  21.40 1.37 1.35  19.14  23.59 1.00      921     1190
##    m1[10]  22.97  22.98 1.62 1.60  20.29  25.54 1.00      736     1035
##    y1[1]    8.84   8.82 4.04 3.85   2.26  15.57 1.00     1165     1562
##    y1[2]   10.43  10.51 4.05 3.89   3.87  17.15 1.00     1254     1649
##    y1[3]   12.02  12.02 3.89 3.80   5.54  18.29 1.00     1800     1835
##    y1[4]   13.59  13.66 3.95 3.73   7.10  20.20 1.00     1856     1858
##    y1[5]   15.08  15.07 3.88 3.75   8.62  21.33 1.00     1899     1763
##    y1[6]   16.60  16.60 3.80 3.62  10.50  22.83 1.00     1790     1750
##    y1[7]   18.26  18.20 4.02 4.02  11.92  24.52 1.00     2078     1825
##    y1[8]   19.72  19.76 3.98 3.87  13.06  25.98 1.00     1986     1657
##    y1[9]   21.46  21.53 4.03 3.91  14.70  28.01 1.00     1893     1856
##    y1[10]  22.88  22.94 4.12 3.89  15.99  29.64 1.00     1827     1651

mdl=cmdstan_model('./ex10.stan') 
fn(mdl,data)
## Initial log joint probability = -154.178 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       17      -18.0071    0.00122729   0.000938315      0.9814      0.9814       27    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -18.01
##    b0         4.28
##    b1         2.22
##    s          0.77
##    m0[1]     10.95
##    m0[2]     14.46
##    m0[3]     14.85
##    m0[4]     21.96
##    m0[5]     16.73
##    m0[6]     19.25
##    m0[7]     13.43
##    m0[8]     19.34
##    m0[9]      8.95
##    m0[10]    12.88
##    m0[11]     9.98
##    m0[12]    12.53
##    m0[13]    24.27
##    m0[14]     6.91
##    m0[15]    23.94
##    m0[16]     9.95
##    m0[17]    23.56
##    m0[18]    18.95
##    m0[19]    17.33
##    m0[20]    14.88
##    y0[1]     10.62
##    y0[2]     15.25
##    y0[3]     15.28
##    y0[4]     21.34
##    y0[5]     16.40
##    y0[6]     19.13
##    y0[7]     13.34
##    y0[8]     19.76
##    y0[9]      6.52
##    y0[10]    12.88
##    y0[11]     9.29
##    y0[12]    14.90
##    y0[13]    24.43
##    y0[14]     8.82
##    y0[15]    24.94
##    y0[16]    10.16
##    y0[17]    23.90
##    y0[18]    18.99
##    y0[19]    18.60
##    y0[20]    14.41
##    m1[1]      6.91
##    m1[2]      8.84
##    m1[3]     10.77
##    m1[4]     12.69
##    m1[5]     14.62
##    m1[6]     16.55
##    m1[7]     18.48
##    m1[8]     20.41
##    m1[9]     22.34
##    m1[10]    24.27
##    y1[1]      6.00
##    y1[2]      8.16
##    y1[3]     11.67
##    y1[4]     17.05
##    y1[5]     14.56
##    y1[6]     16.14
##    y1[7]     19.30
##    y1[8]     19.88
##    y1[9]     23.76
##    y1[10]    23.95
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median      sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -19.86 -19.45    1.40 1.04 -22.71 -18.41 1.00      641      821
##    b0       4.37   4.37    0.81 0.74   3.11   5.71 1.01      364      392
##    b1       2.16   2.18    0.16 0.15   1.89   2.41 1.01      422      629
##    s        0.98   0.94    0.31 0.28   0.57   1.53 1.00      852      748
##    m0[1]   10.86  10.87    0.44 0.39  10.13  11.54 1.01      496      625
##    m0[2]   14.27  14.30    0.36 0.35  13.63  14.81 1.01     1103     1116
##    m0[3]   14.65  14.68    0.36 0.36  14.01  15.19 1.01     1249     1234
##    m0[4]   21.56  21.65    0.65 0.64  20.38  22.49 1.00     1064     1208
##    m0[5]   16.48  16.52    0.40 0.40  15.75  17.06 1.01     1968     1422
##    m0[6]   18.93  19.00    0.50 0.50  18.02  19.66 1.00     1496     1405
##    m0[7]   13.27  13.29    0.37 0.35  12.63  13.84 1.01      816     1030
##    m0[8]   19.01  19.08    0.51 0.50  18.10  19.75 1.00     1472     1384
##    m0[9]    8.91   8.92    0.53 0.47   8.06   9.75 1.01      426      520
##    m0[10]  12.74  12.75    0.38 0.36  12.09  13.33 1.01      703      864
##    m0[11]   9.91   9.92    0.48 0.43   9.12  10.66 1.01      453      547
##    m0[12]  12.40  12.41    0.39 0.36  11.74  13.00 1.01      646      803
##    m0[13]  23.80  23.90    0.80 0.78  22.37  24.92 1.00      871     1133
##    m0[14]   6.93   6.93    0.65 0.58   5.89   7.96 1.01      390      415
##    m0[15]  23.49  23.59    0.78 0.75  22.10  24.58 1.00      899     1162
##    m0[16]   9.88   9.89    0.48 0.44   9.09  10.63 1.01      452      547
##    m0[17]  23.11  23.21    0.75 0.73  21.76  24.17 1.00      936     1162
##    m0[18]  18.63  18.70    0.49 0.49  17.76  19.34 1.00     1579     1354
##    m0[19]  17.06  17.11    0.42 0.42  16.30  17.67 1.00     1974     1423
##    m0[20]  14.68  14.71    0.36 0.35  14.04  15.22 1.01     1265     1234
##    y0[1]   11.06  10.90   22.32 1.51   5.30  17.29 1.00     1777     1948
##    y0[2]   12.12  14.27   59.72 1.50   8.01  19.74 1.00     1801     2059
##    y0[3]   12.71  14.66   51.74 1.53   8.79  21.04 1.00     1995     1822
##    y0[4]   22.08  21.59   31.86 1.63  15.55  28.10 1.00     1785     1957
##    y0[5]   19.41  16.59   87.47 1.49  11.40  23.12 1.00     1872     1785
##    y0[6]   19.16  18.98   31.37 1.60  12.45  25.25 1.00     1956     1851
##    y0[7]   12.99  13.30   13.93 1.52   7.26  19.50 1.00     1845     1887
##    y0[8]   12.29  19.04  291.05 1.60  12.60  25.92 1.00     2082     1664
##    y0[9]    8.69   8.87   59.04 1.59   2.86  14.70 1.00     1754     1839
##    y0[10]  12.12  12.73   25.18 1.51   5.66  18.23 1.00     1905     1891
##    y0[11]  12.87   9.94  197.15 1.58   4.12  16.36 1.00     1891     1678
##    y0[12]  13.57  12.39   76.38 1.46   4.81  18.67 1.00     2024     1794
##    y0[13]  17.83  23.88  256.85 1.77  17.83  30.27 1.00     1738     1676
##    y0[14]   6.22   6.97   52.92 1.69   0.51  14.66 1.00     1593     1892
##    y0[15]  23.11  23.57   22.65 1.70  18.29  29.30 1.00     1832     1880
##    y0[16]  11.25   9.92   93.64 1.61   2.12  16.64 1.00     1757     2015
##    y0[17]  19.23  23.17  166.73 1.77  17.31  29.32 1.00     1707     1885
##    y0[18]  18.78  18.67   14.50 1.47  12.54  24.78 1.00     2094     1959
##    y0[19]  15.70  17.11   67.69 1.50  10.71  23.64 1.00     2102     1982
##    y0[20]  17.56  14.76   47.88 1.49   9.04  21.81 1.00     2029     2041
##    m1[1]    6.93   6.93    0.65 0.58   5.89   7.96 1.01      390      415
##    m1[2]    8.80   8.81    0.54 0.48   7.94   9.65 1.01      424      520
##    m1[3]   10.68  10.69    0.45 0.40   9.93  11.38 1.01      486      614
##    m1[4]   12.55  12.57    0.38 0.36  11.90  13.15 1.01      671      825
##    m1[5]   14.43  14.46    0.36 0.35  13.79  14.97 1.01     1158     1144
##    m1[6]   16.30  16.34    0.39 0.40  15.59  16.88 1.01     1937     1380
##    m1[7]   18.18  18.24    0.47 0.47  17.33  18.83 1.00     1715     1340
##    m1[8]   20.05  20.13    0.56 0.55  19.04  20.85 1.00     1254     1343
##    m1[9]   21.93  22.02    0.68 0.66  20.70  22.87 1.00     1030     1163
##    m1[10]  23.80  23.90    0.80 0.78  22.37  24.92 1.00      871     1133
##    y1[1]    8.44   6.91   50.65 1.58   1.01  12.73 1.00     1778     1858
##    y1[2]    8.41   8.77   27.55 1.54   2.75  14.16 1.00     1672     1894
##    y1[3]   16.95  10.69  274.03 1.47   4.35  16.88 1.00     1935     1750
##    y1[4]  -19.11  12.54 1210.82 1.44   4.63  17.09 1.00     1838     1440
##    y1[5]   14.83  14.46   24.95 1.49   7.85  21.66 1.00     1940     1973
##    y1[6]   13.39  16.31   74.00 1.57   9.63  22.88 1.00     1748     1994
##    y1[7]   18.47  18.25   31.24 1.61  11.43  24.74 1.00     1857     2033
##    y1[8]   22.11  20.10   60.92 1.67  13.71  26.85 1.00     1985     2017
##    y1[9]   21.26  21.96   33.49 1.61  16.26  28.30 1.00     1852     1427
##    y1[10]  24.19  23.84   33.26 1.78  17.41  30.22 1.00     1892     2014

censored data

objective variable has NA

(x,y) i=1-N
(x0,y0) i=1-N0
x1 i=1-N1, y1=NA
(x,y)~N((mx,my),(sx2,sy2,sxy))
(x0,y0)~N((mx,my),(sx2,sy2,sxy))
x1~N(mx,sx2)

ex11-0.stan

data{
  int N0;
  array[N0] vector[2] xy;
  int N1;
  vector[N1] x1;
}
parameters{
  vector[2] m;
  cov_matrix[2] s;
}
model{
  target+=multi_normal_lpdf(xy | m, s);
  x1~normal(m[1],s[1,1]^.5);
}
generated quantities{
  vector[2] xy1;
  xy1=multi_normal_rng(m,s);
  real cr;
  cr=s[1,2]/(s[1,1]*s[2,2])^.5;
}
n=30
x=runif(n,0,9)
y=rnorm(n,10+3*x,4)
cor(x,y)
## [1] 0.899
qplot(x,y)

L=4
n0=sum(x>L)
x0=x[x>L]
y0=y[x>L]
x1=x[x<=L]
n1=sum(x<=L)
cor(x0,y0)
## [1] 0.611
qplot(x0,y0)

mdl=cmdstan_model('./ex11-0.stan') 

data=list(N0=n0,xy=matrix(c(x0,y0),ncol=2),N1=n1,x1=x1)

mle=mdl$optimize(data=data)
## Initial log joint probability = -11806.3 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       70      -102.836    0.00249015    0.00278684           1           1       84    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__    -102.84
##    m[1]       4.48
##    m[2]      19.81
##    s[1,1]     9.14
##    s[2,1]    36.64
##    s[1,2]    36.64
##    s[2,2]   169.72
##    xy1[1]    10.96
##    xy1[2]    42.13
##    cr         0.93
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.3 seconds.
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 finished in 0.5 seconds.
## Chain 3 finished in 0.4 seconds.
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 finished in 0.5 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.4 seconds.
## Total execution time: 0.6 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median     sd    mad      q5    q95 rhat ess_bulk ess_tail
##    lp__   -97.54 -97.20   1.70   1.50 -100.78 -95.51 1.00      430      727
##    m[1]     4.53   4.52   0.61   0.62    3.57   5.55 1.01      590      777
##    m[2]    21.26  21.17   5.92   5.75   11.45  31.29 1.01      136      103
##    s[1,1]  11.75  11.24   3.52   3.20    7.22  18.09 1.00      795     1027
##    s[2,1]  42.02  40.74  24.01  21.39    3.44  80.93 1.02      159      230
##    s[1,2]  42.02  40.74  24.01  21.39    3.44  80.93 1.02      159      230
##    s[2,2] 226.14 190.22 151.62 128.70   51.37 509.45 1.02      144      195
##    xy1[1]   4.58   4.64   3.55   3.49   -1.46  10.27 1.00     1790     1846
##    xy1[2]  21.26  22.47  15.84  13.87   -6.84  44.80 1.00      923      578
##    cr       0.81   0.91   0.27   0.07    0.15   0.97 1.01      199      251

xy=mcmc$draws('xy1')
cor(xy[,,1],xy[,,2])
## [1] 0.764
qplot(xy[,,1],xy[,,2])

objective variable is censored

y i=1-N, y~N(m,s)  
  actual          ya i=1-Na
  lower censored  yl i=1-Nl, y<L, P(y<L)=cdf(L-m /s)
  upper censored  yu i=1-Nu, y>U, P(y>U)=ccdf(U-m /s)

cdf(z) cumulative normal density function, P((-Infinity,z],z~N(0,1))
ccdf(z) complementary CDF, P([z,Infinity),z~N(0,1))

P(y | x,m,s)=P(ya i=1-Na)* P(yl i=1-Nl)* P(yu i=1-Nu)

ex11-1.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  real L;
  int Nl;
  vector[Nl] xl;
  real U;
  int Nu;
  vector[Nu] xu;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
  for(i in 1:Nl)
    target+=normal_lcdf(L | b0+b1*xl[i],s);
  for(i in 1:Nu)
    target+=normal_lccdf(U | b0+b1*xu[i],s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n0=20
x=runif(n0,0,9)
y=rnorm(n0,10+3*x,4)
L=15
y[y<L]=L
nl=sum(y==L)
U=30
y[y>U]=U
nu=sum(y==U)
n=n0-nl-nu
qplot(x,y)

xy0=tibble(x=x,y=y)
xya=filter(xy0, y>L & y<U)
xyl=filter(xy0, y==L)
xyu=filter(xy0, y==U)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
mdl=cmdstan_model('./ex8-0.stan')

data=list(N=n,x=xya$x,y=xya$y,N1=n1,x1=x1)

mle=mdl$optimize(data=data)
## Initial log joint probability = -337.674 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       32      -13.2188   0.000355531   7.26668e-06           1           1       45    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -13.22
##    b0        15.37
##    b1         1.74
##    s          2.63
##    m0[1]     22.77
##    m0[2]     20.67
##    m0[3]     26.59
##    m0[4]     23.85
##    m0[5]     26.79
##    m0[6]     17.83
##    m0[7]     23.75
##    m0[8]     27.31
##    m0[9]     22.03
##    y0[1]     23.68
##    y0[2]     19.18
##    y0[3]     24.92
##    y0[4]     24.52
##    y0[5]     28.71
##    y0[6]     15.77
##    y0[7]     25.93
##    y0[8]     26.79
##    y0[9]     19.70
##    m1[1]     15.74
##    m1[2]     17.40
##    m1[3]     19.05
##    m1[4]     20.70
##    m1[5]     22.35
##    m1[6]     24.01
##    m1[7]     25.66
##    m1[8]     27.31
##    m1[9]     28.97
##    m1[10]    30.62
##    y1[1]     14.98
##    y1[2]     17.36
##    y1[3]     17.70
##    y1[4]     21.34
##    y1[5]     22.11
##    y1[6]     27.50
##    y1[7]     23.26
##    y1[8]     26.99
##    y1[9]     25.34
##    y1[10]    33.60
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,'m')
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -14.23 -13.74 1.76 1.33 -17.83 -12.43 1.01      284      323
##    b0      15.03  15.07 4.58 3.44   8.18  22.24 1.02      264      218
##    b1       1.81   1.80 0.90 0.69   0.38   3.25 1.02      275      246
##    s        3.79   3.40 1.55 1.11   2.22   6.65 1.00      584      676
##    m0[1]   22.71  22.74 1.45 1.19  20.37  24.88 1.01      606      658
##    m0[2]   20.54  20.59 2.11 1.65  17.23  23.90 1.01      319      287
##    m0[3]   26.68  26.62 1.98 1.63  23.66  30.09 1.01      642      740
##    m0[4]   23.84  23.87 1.35 1.15  21.72  25.97 1.01     1280      979
##    m0[5]   26.89  26.81 2.05 1.70  23.76  30.41 1.01      604      654
##    m0[6]   17.58  17.66 3.38 2.60  12.45  22.92 1.02      272      220
##    m0[7]   23.73  23.76 1.35 1.15  21.61  25.83 1.01     1205      950
##    m0[8]   27.43  27.33 2.26 1.88  24.04  31.28 1.01      529      575
##    m0[9]   21.95  21.99 1.63 1.32  19.33  24.47 1.01      430      389
##    y0[1]   22.73  22.69 4.26 3.71  15.96  29.39 1.00     1839     1517
##    y0[2]   20.59  20.79 4.62 3.97  13.23  27.60 1.00      925      838
##    y0[3]   26.75  26.82 4.39 3.75  19.65  33.76 1.00     1369     1430
##    y0[4]   23.75  23.81 4.37 3.63  16.76  30.47 1.00     1828     1773
##    y0[5]   26.81  26.82 4.58 3.84  19.41  34.00 1.00     1241     1154
##    y0[6]   17.56  17.49 5.05 4.23   9.44  25.78 1.00      562      756
##    y0[7]   23.61  23.58 4.34 3.58  16.75  30.66 1.00     1633     1394
##    y0[8]   27.38  27.45 4.59 4.11  19.93  34.44 1.00     1418     1393
##    y0[9]   21.94  21.87 4.36 3.70  15.00  28.81 1.00     1475     1328
##    m1[1]   15.42  15.46 4.39 3.34   8.89  22.27 1.02      265      223
##    m1[2]   17.14  17.18 3.59 2.74  11.75  22.77 1.02      270      219
##    m1[3]   18.85  18.91 2.81 2.15  14.54  23.34 1.02      282      235
##    m1[4]   20.57  20.62 2.10 1.64  17.28  23.91 1.01      320      294
##    m1[5]   22.28  22.32 1.54 1.26  19.78  24.67 1.01      489      491
##    m1[6]   24.00  24.02 1.35 1.16  21.86  26.18 1.00     1393     1055
##    m1[7]   25.72  25.70 1.66 1.38  23.08  28.53 1.01      966      809
##    m1[8]   27.43  27.33 2.27 1.88  24.04  31.28 1.01      528      575
##    m1[9]   29.15  29.02 3.00 2.49  24.60  34.16 1.01      409      451
##    m1[10]  30.86  30.69 3.78 3.06  25.09  37.15 1.01      359      380
##    y1[1]   15.39  15.37 5.94 4.74   6.43  24.41 1.01      431      491
##    y1[2]   17.20  17.32 5.47 4.48   8.37  26.16 1.01      490      564
##    y1[3]   18.90  18.93 4.88 4.10  11.40  26.59 1.00      695      758
##    y1[4]   20.54  20.71 4.62 3.79  13.15  27.37 1.00     1068     1004
##    y1[5]   22.37  22.48 4.32 3.63  15.38  29.01 1.00     1586     1601
##    y1[6]   23.88  23.88 4.38 3.64  17.23  30.76 1.00     2033     1817
##    y1[7]   25.59  25.67 4.36 3.50  18.77  32.36 1.00     1969     1722
##    y1[8]   27.60  27.53 4.83 4.04  20.09  35.41 1.00     1270     1108
##    y1[9]   29.28  29.13 4.89 4.13  21.78  37.01 1.00      984     1277
##    y1[10]  30.78  30.82 5.63 4.62  22.27  39.42 1.01      641      798

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

data=list(N=n,x=xya$x,y=xya$y,
          L=L,Nl=nl,xl=xyl$x,
          U=U,Nu=nu,xu=xyu$x,
          N1=n1,x1=x1)
mdl=cmdstan_model('./ex11-1.stan') 

mle=mdl$optimize(data=data)
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Initial log joint probability = -166.149 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
##       24      -22.7103    0.00150634   0.000243648      0.9838      0.9838       33    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -22.71
##    b0         8.96
##    b1         3.28
##    s          4.43
##    m0[1]     22.89
##    m0[2]     18.95
##    m0[3]     30.08
##    m0[4]     24.93
##    m0[5]     30.46
##    m0[6]     13.59
##    m0[7]     24.74
##    m0[8]     31.44
##    m0[9]     21.50
##    y0[1]     25.16
##    y0[2]     17.33
##    y0[3]     36.80
##    y0[4]     19.44
##    y0[5]     30.77
##    y0[6]     15.86
##    y0[7]     30.83
##    y0[8]     32.83
##    y0[9]     21.00
##    m1[1]      9.67
##    m1[2]     12.78
##    m1[3]     15.89
##    m1[4]     19.00
##    m1[5]     22.11
##    m1[6]     25.22
##    m1[7]     28.33
##    m1[8]     31.45
##    m1[9]     34.56
##    m1[10]    37.67
##    y1[1]     15.56
##    y1[2]     21.52
##    y1[3]     16.84
##    y1[4]     18.80
##    y1[5]     21.75
##    y1[6]     22.76
##    y1[7]     29.49
##    y1[8]     30.55
##    y1[9]     29.59
##    y1[10]    28.64
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,'m',ch=T)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -23.11 -22.71 1.55 1.23 -26.02 -21.44 1.02      287      481
##    b0       5.84   6.63 5.40 4.43  -3.62  12.30 1.02      220      281
##    b1       3.99   3.80 1.14 0.93   2.63   5.92 1.01      240      255
##    s        6.51   5.94 2.66 1.95   3.71  11.17 1.01      405      602
##    m0[1]   22.81  22.91 1.91 1.70  19.64  25.65 1.00      751      670
##    m0[2]   18.01  18.36 2.48 2.11  13.64  21.21 1.01      314      413
##    m0[3]   31.57  31.15 2.97 2.50  27.72  36.81 1.01      641      526
##    m0[4]   25.29  25.24 1.96 1.76  22.37  28.45 1.00     1344      882
##    m0[5]   32.03  31.61 3.07 2.58  28.06  37.49 1.01      604      504
##    m0[6]   11.48  12.07 3.94 3.27   4.41  16.17 1.01      229      287
##    m0[7]   25.06  25.02 1.95 1.76  22.14  28.15 1.00     1290      938
##    m0[8]   33.22  32.72 3.34 2.79  28.97  39.27 1.01      529      490
##    m0[9]   21.11  21.30 2.03 1.77  17.67  23.96 1.01      508      523
##    y0[1]   23.03  22.93 7.07 5.98  12.04  34.19 1.00     1770     1714
##    y0[2]   17.97  18.27 7.47 6.37   5.61  28.79 1.00     1334      973
##    y0[3]   31.35  30.98 7.66 6.55  19.90  43.92 1.01     1782     1284
##    y0[4]   24.99  24.94 7.09 6.20  14.32  36.51 1.00     1983     1576
##    y0[5]   32.06  31.63 7.45 6.10  20.99  44.27 1.00     1701     1715
##    y0[6]   11.19  12.10 8.06 6.68  -2.74  21.93 1.00      658      965
##    y0[7]   25.03  25.01 7.28 6.18  13.60  36.42 1.00     2005     1688
##    y0[8]   33.22  32.78 7.88 6.66  21.63  46.19 1.00     1605      968
##    y0[9]   21.22  21.41 7.37 6.01   9.65  32.51 1.00     1833     1680
##    m1[1]    6.71   7.48 5.17 4.25  -2.33  12.86 1.02      221      276
##    m1[2]   10.50  11.12 4.18 3.49   3.12  15.44 1.02      226      294
##    m1[3]   14.29  14.79 3.26 2.72   8.32  18.29 1.01      244      324
##    m1[4]   18.07  18.42 2.46 2.11  13.73  21.27 1.01      316      418
##    m1[5]   21.86  22.02 1.96 1.71  18.57  24.69 1.01      598      564
##    m1[6]   25.65  25.59 1.99 1.79  22.70  28.88 1.00     1413      857
##    m1[7]   29.44  29.12 2.53 2.23  26.05  33.79 1.01      940      733
##    m1[8]   33.23  32.72 3.34 2.79  28.98  39.28 1.01      528      490
##    m1[9]   37.02  36.29 4.28 3.55  31.68  44.70 1.01      404      450
##    m1[10]  40.81  39.91 5.26 4.25  34.31  50.17 1.01      351      397
##    y1[1]    6.82   7.81 8.84 7.44  -8.51  18.62 1.01      596      590
##    y1[2]   10.33  10.85 8.36 6.66  -2.84  21.60 1.01      692      685
##    y1[3]   14.57  14.98 7.78 6.79   1.78  26.13 1.00     1172     1116
##    y1[4]   18.06  18.30 7.67 6.07   6.18  29.56 1.00     1775     1200
##    y1[5]   22.03  22.35 7.32 6.40  10.10  33.31 1.00     1714     1227
##    y1[6]   25.58  25.49 7.19 6.04  14.84  37.09 1.00     1947     1693
##    y1[7]   29.23  29.03 7.32 6.27  18.33  41.06 1.00     2056     1366
##    y1[8]   33.18  32.64 7.47 6.40  21.99  46.14 1.00     1692     1351
##    y1[9]   37.24  36.42 8.56 6.99  24.87  51.54 1.01     1213     1103
##    y1[10]  40.54  39.77 8.42 7.04  28.84  54.86 1.00      785      988

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

sensitivity/specificity analysis

sensitivity: true positive rate TPR = TP/(TP+FN)
specificity: true negative rate TFR = TN/(FP+TN)
ROC curve: se vs 1-sp

positive predictive value ppv = TP/(TP+FP)
negative predictive value npv = TN/(TN+FN)

ex14.stan

estimating sens and spec

data {
  int N;
  array[N] int x;
  array[N] int y;
}
parameters {
  real<lower=0,upper=1> p;
  real<lower=0,upper=1> se;
  real<lower=0,upper=1> sp;
}
model {
  p~uniform(0,1);
  se~uniform(0,1);
  sp~uniform(0,1);
  for (i in 1:N) {
    y[i]~bernoulli(x[i]*se+(1-x[i])*(1-sp));
  }
}
generated quantities {
  real ppv;
  real npv;
  ppv=se*p/((se*p)+((1-p)*(1-sp)));
  npv=(1-p)*sp/(((1-p)*sp)+(p*(1-se)));
}
n=20
x=sample(0:1,n,replace=T)
p=(x+0.5)*0.5
y=rbinom(n,1,p)
data=list(N=n,x=x,y=y)

mdl=cmdstan_model('./ex14.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -14.7148 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        6      -13.4428   0.000615617   0.000249883           1           1        9    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -13.44
##      p        0.88
##      se       0.58
##      sp       0.63
##      ppv      0.92
##      npv      0.17
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -19.30 -18.95 1.32 1.03 -21.77 -17.89 1.00      755      840
##      p      0.49   0.48 0.28 0.36   0.06   0.94 1.01      710      423
##      se     0.57   0.58 0.13 0.14   0.35   0.77 1.00     3177     1523
##      sp     0.59   0.61 0.15 0.16   0.34   0.82 1.00     1820     1245
##      ppv    0.55   0.58 0.29 0.37   0.08   0.96 1.01      742      556
##      npv    0.57   0.60 0.28 0.36   0.07   0.96 1.01      749      594

ppv=mcmc$draws('ppv')
npv=mcmc$draws('npv')

qplot(ppv,npv)

2x2 cross table

Effect occur y=1 with probabilty p01, p11 from each cause x{0,1}
event frequncy nxy of effect y{0,1} by cause x{0,1}
n01~B(n0.,p0)
n11~B(n1.,p1)

n01=n0p0, n00=n0(1-p0)
n11=n1p1, n10=n1(1-p1)

p00=n00/n=n0(1-p0)/n, p01=n01/n=n0p0/n
p10=n10/n=n1(1-p1)/n, p11=n11/n=n1p1/n

Cramer'V  (chi2/n/(min(row,column)-1))^.5
  in 2x2  
  crv =(n11*n00-n10*n01)/(n0.*n1.*n.0*n.1)^.5
      =(n0(1-p0)n1p1-n0p0n1(1-p1))/(n0n1(n0(1-p0)+n1(1-p1))(n0p0n1p1))^.5

kappa coefficient   k=(po-pe)/(1-pe)
  po: Observed agreement (proportion of times both raters agreed)
  pe: Expected agreement under independence
      po=p00+p11
        =(n0(1-p0)+n1p1)/n
      pe=(p00+p01)(p00+p10)(p10+p11)(p01+p11)
        =n0/n*(n0(1-p0)+n1(1-p1))/n*(n0p0+n1p1)/n*n1/n

ex16-1.stan

data {
  int n;
  int n0;
  int n01;
  int n1;
  int n11;
}
parameters {
  real<lower=0,upper=1> p0;
  real<lower=0,upper=1> p1;
}
model {
  n01~binomial(n0,p0);
  n11~binomial(n1,p1);
}
generated quantities {
  real RR;
  RR=p1/p0;
  real OR;
  OR=(p1/(1-p1))/(p0/(1-p0));
}

ex16-2.stan

data {
  int n;
  int n0;
  int n01;
  int n1;
  int n11;
}
parameters {
  real<lower=0,upper=1> p0;
  real<lower=0,upper=1> p1;
}
model {
  n01~binomial(n0,p0);
  n11~binomial(n1,p1);
}
generated quantities {
  real RR;
  RR=p1/p0;
  real OR;
  OR=(p1/(1-p1))/(p0/(1-p0));
  real crv;
  crv=(n0*(1-p0)*n1*p1-n0*p0*n1*(1-p1))/(n0*n1*(n0*(1-p0)+n1*(1-p1))*(n0*p0+n1*p1))^.5;
  real k;
  real po;
  po=(n0*(1-p0)+n1*p1)/n;
  real pe;
  pe=n0/n*(n0*(1-p0)+n1*(1-p1))/n*(n0*p0+n1*p1)/n*n1/n;
  k=(po-pe)/(1-pe);
}
n0=30
n01=rbinom(1,n0,0.3)
n1=30
n11=rbinom(1,n1,0.6)
data=list(n=n0+n1,n0=n0,n01=n01,n1=n1,n11=n11)

mdl=cmdstan_model('./ex16-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -55.3333 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        5      -40.7173    0.00355404   6.49073e-05           1           1        8    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -40.72
##      p0       0.40
##      p1       0.43
##      RR       1.08
##      OR       1.15
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -44.48 -44.20 0.89 0.66 -46.25 -43.58 1.00      939      972
##      p0     0.40   0.40 0.08 0.08   0.27   0.54 1.00     2047     1583
##      p1     0.44   0.44 0.08 0.09   0.30   0.58 1.00     2165     1521
##      RR     1.13   1.08 0.33 0.31   0.67   1.73 1.00     2068     1567
##      OR     1.30   1.14 0.67 0.55   0.51   2.52 1.00     2109     1556

data=list(n=n0+n1,n0=n0,n01=n01,n1=n1,n11=n11)

mdl=cmdstan_model('./ex16-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -58.2846 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        4      -40.7173     0.0154814   0.000587882      0.9942      0.9942        7    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -40.72
##      p0       0.40
##      p1       0.43
##      RR       1.08
##      OR       1.15
##      crv      0.03
##      k        0.49
##      po       0.52
##      pe       0.06
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -44.48 -44.20 0.89 0.66 -46.25 -43.58 1.00      939      972
##      p0     0.40   0.40 0.08 0.08   0.27   0.54 1.00     2047     1583
##      p1     0.44   0.44 0.08 0.09   0.30   0.58 1.00     2165     1521
##      RR     1.13   1.08 0.33 0.31   0.67   1.73 1.00     2068     1567
##      OR     1.30   1.14 0.67 0.55   0.51   2.52 1.00     2109     1556
##      crv    0.03   0.03 0.12 0.12  -0.17   0.22 1.00     2111     1504
##      k      0.49   0.48 0.06 0.06   0.38   0.58 1.00     2111     1553
##      po     0.52   0.52 0.06 0.06   0.42   0.61 1.00     2115     1628
##      pe     0.06   0.06 0.00 0.00   0.06   0.06 1.00     2065     1609

point of subjective equality PSE

PSE: 50% threshold for sensing the difference between two stimuli is equal
JND: Just noticeable difference, difference between 50% threshold and 75%

r~B(n,p) #reaction for stimuli
p=1/(1+exp(-(a+b*x)))
x=x1-x0, x0,x1 is strength of stimuli

PSE=-a/b
JND=(log(0.75/0.25)-a)/b-PSE

ex6-3-0.stan

mulit logistic regression

data{
  int N;
  int m;
  vector[N] x;
  array[N] int y;
}
parameters{
  real b0;
  real b1;
}
model{
  y~binomial_logit(m,b0+b1*x);
}
n=20
m=10
x=runif(n,-2,2)
y=rbinom(n,m,1/(1+exp(-(-2+3*x))))

glm(y/m~x,family=binomial('logit'))
## 
## Call:  glm(formula = y/m ~ x, family = binomial("logit"))
## 
## Coefficients:
## (Intercept)            x  
##       -2.32         2.78  
## 
## Degrees of Freedom: 19 Total (i.e. Null);  18 Residual
## Null Deviance:       16.6 
## Residual Deviance: 1.9   AIC: 10
data=list(N=n,m=m,x=x,y=y)

mdl=cmdstan_model('./ex6-3-0.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -286.995 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       10      -52.0684   0.000462636    0.00031605           1           1       12    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -52.07
##      b0      -2.32
##      b1       2.78
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -53.07 -52.77 0.98 0.71 -55.09 -52.12 1.00      751      791
##      b0    -2.44  -2.42 0.44 0.46  -3.22  -1.76 1.00      670      688
##      b1     2.92   2.88 0.43 0.43   2.29   3.71 1.00      660      722

b0=mcmc$draws('b0') |> as.vector()
b1=mcmc$draws('b1') |> as.vector()

pse=-b0/b1
quantile(pse,probs=c(0.0,0.025,0.05,0.25,0.5,0.75,0.95,0.975,1),na.rm=T)
##    0%  2.5%    5%   25%   50%   75%   95% 97.5%  100% 
## 0.531 0.665 0.690 0.778 0.837 0.896 0.981 1.009 1.196
jnd=(log(0.75/0.25)-b0)/b1-pse
quantile(jnd,probs=c(0.0,0.025,0.05,0.25,0.5,0.75,0.95,0.975,1),na.rm=T)
##    0%  2.5%    5%   25%   50%   75%   95% 97.5%  100% 
## 0.247 0.287 0.296 0.344 0.382 0.420 0.480 0.508 0.694
x1=runif(length(b0),-2,2)
p=1/(1+exp(-(b0+b1*x1)))
pm=1/(1+exp(-(median(b0)+median(b1)*x1)))

qplot(x1,pm,col=I('darkgray'),ylab='p')+
  geom_line(aes(x=x1,p=p),col='red')